{
 "metadata": {
  "language": "Julia",
  "name": "",
  "signature": "sha256:00e0fe30cd2c454878301c4e1accd538c844a6438a23c55e0d332ecebb33136f"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Polyphase Filtering Explained"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Glossery ##\n",
      "\n",
      "The folowing list contains naming conventions used in this document and Mulirate's source code. You will see references to columns and rows. Take note that Julia is column major language.\n",
      "\n",
      "* `h`: a vector of filter coefficients (aka taps)\n",
      "* `PFB`: poly phase filter bank\n",
      "* `<symbol>Len`: a quantity of something\n",
      "* `\ud835\udf19`: a phase of a `PFB`\n",
      "* `\u0192`: frequency\n",
      "* `<symbol>Time`: vector of time values associated with another vector <symbol>\n",
      "\n",
      "## Background ##\n",
      "\n",
      "A conventional polyphase FIR filter implementation takes filter taps and splits them up into phases. Why are they called phases? I still don't know. But the important thing is that this allows for efficient interpolation. In the naive approach to interpolation by a factor of `L` one inserts, or 'stuffs', `L-1` zeros between each input sample prior to lowpass filtering:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "interpolation = 4\n",
      "xLen          = 4\n",
      "x             = [1:xLen]\n",
      "xStuffed      = vec(vcat( zeros(Int,3,4), x' ))\n",
      "println( \"x = $(x')\" )\n",
      "println( \"xStuffed = $(xStuffed)\" )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "x = [1 2 3 4]\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "xStuffed = [0,0,0,1,0,0,0,2,0,0,0,3,0,0,0,4]\n"
       ]
      }
     ],
     "prompt_number": 78
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This inefficient however. The zeros that we stuffed inbetween the input samples contribute nothing to the output of the filtering process other than placeholders to keep filter taps properly aligned with input. Remeber that FIR filter involves convolving `x` and `h`. Let's compute, using strings for readablity, the first four output samples with hypothetical vector of taps `h`:\n",
      "\n",
      "    y[1] = x1*h12 + x2*h11 + x3*h10 + x4*h9 + x5*h8 + x6*h7 + x7*h6 + x8*h5 + x9*h4 + x10*h3 + x11*h2 + x12*h1\n",
      "    y[2] = x2*h12 + x3*h11 + x4*h10 + x5*h9 + x6*h8 + x7*h7 + x8*h6 + x9*h5 + x10*h4 + x11*h3 + x12*h2 + x13*h1\n",
      "    y[3] = x3*h12 + x4*h11 + x5*h10 + x6*h9 + x7*h8 + x8*h7 + x9*h6 + x10*h5 + x11*h4 + x12*h3 + x13*h2 + x14*h1\n",
      "    y[4] = x4*h12 + x5*h11 + x6*h10 + x7*h9 + x8*h8 + x9*h7 + x10*h6 + x11*h5 + x12*h4 + x13*h3 + x14*h2 + x15*h1"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hLen = 12\n",
      "\n",
      "for yIdx in 1:4\n",
      "    hIdx = hLen\n",
      "    print( \"y[$yIdx] = \" )\n",
      "    xStartIdx = yIdx\n",
      "    xStopIdx  = xStartIdx+hLen-1\n",
      "    for xIdx in xStartIdx:xStopIdx\n",
      "#         print( \"$(xStuffed[xIdx])*h$hIdx\" )\n",
      "        print( \"x$xIdx*h$hIdx\" )        \n",
      "        xIdx < xStopIdx && print( \" + \" )\n",
      "        hIdx -= 1\n",
      "    end\n",
      "    println()    \n",
      "end\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "y[1] = x1*h12"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " + x2*h11 + x3*h10 + x4*h9 + x5*h8 + x6*h7 + x7*h6 + x8*h5 + x9*h4 + x10*h3 + x11*h2 + x12*h1\n",
        "y[2] = x2*h12 + x3*h11 + x4*h10 + x5*h9 + x6*h8 + x7*h7 + x8*h6 + x9*h5 + x10*h4 + x11*h3 + x12*h2 + x13*h1\n",
        "y[3] = x3*h12 + x4*h11 + x5*h10 + x6*h9 + x7*h8 + x8*h7 + x9*h6 + x10*h5 + x11*h4 + x12*h3 + x13*h2 + x14*h1\n",
        "y[4] = x4*h12 + x5*h11 + x6*h10 + x7*h9 + x8*h8 + x9*h7 + x10*h6 + x11*h5 + x12*h4 + x13*h3 + x14*h2 + x15*h1\n"
       ]
      }
     ],
     "prompt_number": 79
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See a pettern? Let's make it more obvious by skipping miltiplications where `x == 0`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for yIdx in 1:4\n",
      "    xStartIdx = yIdx\n",
      "    xStopIdx  = xStartIdx+hLen-1\n",
      "    hIdx      = hLen\n",
      "    xPrevZero = true\n",
      "    for xIdx in xStartIdx:xStopIdx\n",
      "        if xStuffed[xIdx] != 0\n",
      "            xPrevZero == false && print( \" + \" )\n",
      "            print( \"$(xStuffed[xIdx])*h$hIdx\" )\n",
      "            xPrevZero = false\n",
      "        end\n",
      "        hIdx -= 1\n",
      "    end\n",
      "    println()    \n",
      "end\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1*h9 + "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2*h5 + 3*h1\n",
        "1*h10 + 2*h6 + 3*h2\n",
        "1*h11 + 2*h7 + 3*h3\n",
        "1*h12 + 2*h8 + 3*h4\n"
       ]
      }
     ],
     "prompt_number": 80
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now the pattern is more obious. For the first four `y`'s we can decompose `h` into four phases, `\ud835\udf191 = [ h1, h5, h9 ]`, `\ud835\udf192 = [ h2, h6, h10 ]`, `\ud835\udf193 = [ h4, h7, h11 ]`, and `\ud835\udf194 = [ h4, h8, h12 ]`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "using Multirate\n",
      "using DSP\n",
      "using PyPlot"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "INFO: Loading help data...\n"
       ]
      }
     ],
     "prompt_number": 81
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Design Finite Impulse Response (FIR) Filter Taps\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "N\ud835\udf19              = 32                                           # Number of polyphase partitions\n",
      "tapsPer\ud835\udf19        = 10                                           # N\ud835\udf19 * tapsPer\ud835\udf19 will be the length of out prototype filter taps\n",
      "resampleRate    = float64(pi)                                  # Can be any arbitrary resampling rate\n",
      "polyorder       = 4                                            # Our taps will tranformed into\n",
      "Nx              = 40                                           # Number of signal samples\n",
      "t               = 0:Nx-1                                       # Time range\n",
      "x\u01921             = 0.15                                         # First singal frequency\n",
      "x\u01922             = 0.3                                          # Second signal frequency\n",
      "x               = cos(2*pi*x\u01921*t) .+ 0.5sin(2*pi*x\u01922*t*pi)     # Create the two signals and add them\n",
      "tx              = [0:length(x)-1]                              # Signal time vector\n",
      "cutoffFreq      = min( 0.45/N\ud835\udf19, resampleRate/N\ud835\udf19 )              # N\ud835\udf19 is also the integer interpolation, so set cutoff frequency accordingly\n",
      "hLen            = tapsPer\ud835\udf19*N\ud835\udf19                                  # Total number of filter taps\n",
      "h               = firdes( hLen, cutoffFreq, DSP.kaiser ) .* N\ud835\udf19 # Generate filter taps and scale by polyphase interpolation (N\ud835\udf19)\n",
      "myfilter        = FIRFilter( h, resampleRate, N\ud835\udf19, polyorder )  # Construct a FIRFilter{FIRFarrow} object\n",
      "y               = filt( myfilter, x )                          # Filter x\n",
      "ty              = [0:length(y)-1]./resampleRate - tapsPer\ud835\udf19/2   # Create y time vector. Accout for filter delay so the plots line up\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 82,
       "text": [
        "126-element Array{Float64,1}:\n",
        " -5.0    \n",
        " -4.68169\n",
        " -4.36338\n",
        " -4.04507\n",
        " -3.72676\n",
        " -3.40845\n",
        " -3.09014\n",
        " -2.77183\n",
        " -2.45352\n",
        " -2.13521\n",
        " -1.8169 \n",
        " -1.49859\n",
        " -1.18028\n",
        "  \u22ee      \n",
        " 31.2873 \n",
        " 31.6056 \n",
        " 31.9239 \n",
        " 32.2423 \n",
        " 32.5606 \n",
        " 32.8789 \n",
        " 33.1972 \n",
        " 33.5155 \n",
        " 33.8338 \n",
        " 34.1521 \n",
        " 34.4704 \n",
        " 34.7887 "
       ]
      }
     ],
     "prompt_number": 82
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}